function [tsOut,IT] = detrendZeroMean(ts,NC)
%detrendZeroMean - Computation of detrended, zero mean time series
%
%       [tsOut,IT] = detrendZeroMean(ts,NC)
% 
% This function computes the detrended, zero mean version of the time
% series related to the tsData object ts for the number of components
% NC (allowed values: 2 for East-North, 3 for East-North-Vertical).
% If the detrending is successful, the output variable tsOut is upgraded
% with the obtained detrendedZeroMeanE,..., properties and IT is true.
% If the detrending cannot be carried out because the trend is not still
% availabe, tsOut = ts and It is false.
% If NC is undefined, empty or is neither 2 nor 3, NC = 3 is used.
%
% See also CMEstackFiltering.

% G. Teza, 2022.

if (nargin < 2) || isempty(NC) || (numel(NC) > 1) || ~ismember(NC,[2,3])
    NC = 3;
end

if NC == 2
    VC = [0 1];
else
    VC = [0 1 2];
end

tsOut = ts;
for h = 1:NC
    if VC(h) == 0
        if ~isempty(ts.estimatedTrendE) 
            t = ts.estimatedTrendE.t;
            o = ts.estimatedTrendE.o;   % data in m
            trend = ts.estimatedTrendE.trend.trend;
            offs = ts.offsetsE;
        else
            IT = false;
            return
        end
    end
    if VC(h) == 1
        if ~isempty(ts.estimatedTrendN)
            t = ts.estimatedTrendN.t;
            o = ts.estimatedTrendN.o;   % data in m
            trend = ts.estimatedTrendN.trend.trend;
            offs = ts.offsetsN;
        end
    end
    if VC(h) == 2
        if ~isempty(ts.estimatedTrendV)
            t = ts.estimatedTrendV.t;
            o = ts.estimatedTrendV.o;   % data in m
            trend = ts.estimatedTrendV.trend.trend;
            offs = ts.offsetsV;
        else
            fprintf('WARNING: station %s, vertical trend unavailable \n',...
                ts.statName);
            continue
        end
    end
    lts = numel(t);
    tf = serial2frac(t);    % to have t in fractional year
    v = trend/1000;         % velocity in m/y
    et = v*(tf-tf(1));      % raw trend
    od = o-et;              % raw residual
    if isempty(offs)
        mod = mean(od);
        odm = od-mod;       % zero mean residual time series
    else
        nos = numel(offs);
        odm = zeros(lts,1);
        for k = 1:(nos+1)
            if k == 1
                t1k = t(1);
            else
                [~,I1k] = min(abs(t-offs(k-1)));
                if t(I1k) <= offs(k-1)  % the offset could be not in t
                    t1k = t(I1k+1);
                else
                    t1k = t(I1k);
                end
            end
            if k == nos+1
                t2k = t(end);
            else
                [~,I2k] = min(abs(t-offs(k)));
                if t(I2k) < offs(k)  % the offset could be not in t
                    t2k = t(I2k);
                else
                    t2k = t(I2k-1);
                end
            end
            dtk = (t1k:t2k)';
            Itk = ismember(t,dtk);  % true for the t in interval dtk
            odk = od(Itk);
            modk = mean(odk);
            odmk = odk-modk;        % detrended, zero-mean in [t1k t2k]
            odm(Itk) = odmk;
        end
    end
    if VC(h) == 0
        tsOut.detrendedZeroMeanE.t = t;
        tsOut.detrendedZeroMeanE.d = odm;
        tsOut.detrendedZeroMeanE.method = 'Modeled';
    end
    if VC(h) == 1
        tsOut.detrendedZeroMeanN.t = t;
        tsOut.detrendedZeroMeanN.d = odm;
        tsOut.detrendedZeroMeanN.method = 'Modeled';
    end
    if VC(h) == 2
        tsOut.detrendedZeroMeanV.t = t;
        tsOut.detrendedZeroMeanV.d = odm;
        tsOut.detrendedZeroMeanN.method = 'Modeled';
    end
end
IT = true;